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Abstract.  The  vorticity- vector  potential  formulation  and  the  fourth- 
order  difference  method  are  combined  with  LES  with  dynamic  subgrid  scale 
model.  The  new  sharp  cut  off  filtering  is  devised,  which  has  the  fourth  order 
accuracy.  The  present  computational  method  is  validated  for  the  simula- 
tion of  transition  in  channel.  The  channel  flow  with  two  rows  of  transverse 
riblet  is  computed  by  the  present  method  and  is  shown  that  the  drag  reduc- 
tion is  attained  after  the  turbulet  field  developes  fully  on  the  riblet  surface. 
It  is  also  found  that  the  transverse  riblet  has  the  same  effect  as  the  suction. 


1.  Introduction 

Recently,  the  large  eddy  simulation  is  devoted  to  study  the  turbulence 
phenomena  with  complex  geometric  configulations  such  as  turbulence  over 
the  airfoil  (Weber  et  al  2000)  and  the  turbulence  flow  in  motored  piston- 
cylinder  assembly  (Verzicco  et  al  2000).  In  conducting  these  simulations 
accurately,  it  is  crucially  important  to  use  the  accurate  large  eddy  simula- 
tion technique  in  the  generalized  coordinate  system.  The  accurate  filtering 
technique  is  also  important  as  well  as  accurate  computational  method.  The 
dynamic  subgrid  scale  model  (Germano  et  al  1991)  should  be  adopted  in 
these  simulations.  It  has  also  to  be  considered  the  importance  of  the  formu- 
lation method  in  dealing  with  incompressible  turbulent  flow.  The  fractional 
step  method  has  been  widely  used  as  the  formulation  method  since  the  reg- 
ular grid  is  applicable  and  the  continuity  equation  of  the  velocity  field  is 
satisfied  sufficiently. 
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Although,  the  vorticity-vector  potential  method  has  been  used  in  the 
early  stage  of  the  large  eddy  simulations  (Mansour  et.  al. , 1979),  this 
method  was  replaced  to  the  fractional  step  method  due  to  the  difficulty 
on  applicating  the  wall  boundary  conditions  adequately.  However,  recently 
this  problem  has  been  dramatically  solved  by  the  research  works  of  the 
present  author  (Tokunaga,  1992;  Shimada  et.  al.,  1991;  Tokunaga  et.  al, 
1991;  Tokunaga,  1999).  The  fourth  order  accurate  difference  method  has 
been  also  combined  to  this  formulation  in  the  frame  work  of  the  dynamic 
subgrid  scale  model. 

The  present  author  conducts  the  research  works  on  the  drag  reduction  of 
the  transverse  riblet  by  making  use  of  the  large  eddy  simulation  technique 
which  consists  of  the  fourth  order  difference  method,  fourth  order  accurate 
sharp-cut-off  filter,  dynamic  subgrid  scale  model  and  the  vorticity-vector 
potential  formaulation  in  the  curved  coordinate  system. 

It  has  to  be  pointed  out  at  first  the  direct  numerical  simulation  over 
riblet  (Choi  et.  al.,  1994)  about  research  works  of  turbulent  drag  reduction. 
However,  another  studies  exist  on  numerical  simulations  of  the  turbulent 
drag  reduction,  such  as  the  linear  feedback  control  (Cortelezzi  et.  al.,  1999) 
as  well  as  the  transverse  riblet  (Tokunaga  and  Yamauchi,  1997;  Tokunaga, 
1999;  Manuilovich,  1994;  Liu  et.  al,  1994;  Mochizuki  et.  al,  1994).  The 
present  author  showed  that  the  concave  transverse  riblet  has  the  signifi- 
cant effect  for  drag  reduction  since  the  streak  structure  of  wall  bounded 
turbulence  is  blocked  by  the  concave  transverse  riblet  and  the  turbulence 
is  weakened  as  the  suction  (Tokunaga,  1999).  The  suction  is  used  as  an 
effective  means  for  laminarizing  the  flow  along  airfoil  and  reducing  the  drag, 
in  which  the  slotted  surface  is  accompanied  for  delaying  the  transition  and 
reducing  the  drag,  although  this  apparatus  includes  the  complex  mechanical 
arrangements  (Bobbit  et.  al,  1991). 

In  the  present  paper,  we  intend  to  verify  more  precisely  the  drag  reduc- 
tion effect  by  conducting  the  large  eddy  simulation  of  the  turbulence  over 
the  two  rows  of  transverse  riblet  in  the  periodic  domain  of  the  streamwise 
direction.  The  computaional  grid  cinsists  of  129  X 85  X 65  where  the  curved 
coordinate  is  used. 

2.  Governing  Equations 

We  adopt  the  fourth  order  accurate  sharp-cut-off  filter  which  is  able  to  be 
applied  to  the  generalized  body  fitted  coordinate  system  since  the  constant 
width  Amin  / 2 is  used  as  shown  in  Fig.  1,  where  Smin/2  denotes  the  minimum 
grid  spacing  in  the  generalized  coordinate,  and  therefore  the  filtering  and 
differentiation  process  is  comutable  each  other  with  respect  to  the  order. 
The  dynamic  subgrid  scale  is  used  and  the  integration  in  the  box  filtering 
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is  performed  analytically  by  approximating  the  turbulent  field  with  the 
second  order  Lagrage  interpolation  both  on  the  basic  and  test  filtering. 
The  integration  is  performed  by  the  Simpson’s  rule,  so  that  the  fourth 
order  accuracy  is  assured. 


Figure  /.  Filtering  process 


The  vorticty-vector  potential  formulation  is  applied  to  the  equation  of 
LES  with  dynamic  subgrid  scale  model.  , The  vorticity  v is  defined  as 

w = v x u.  (1) 

The  vorticity  transport  equation  is  derived  by  taking  the  rotation  of  the 
Navier-Stokes  equation  with  the  SGS  model.  After  introducing  the  vector 
potential 

u = vx^,  (2) 

the  Poisson  equations  are  also  derived  as 

Aipi  = -us,,  (3) 


w 


here  the  solenoidal  condition  is  assumed  for  the  vector  potential 


dj’j 

dxj 


= 0. 


(4) 


These  equations  are  transformed  to  the  generalized  coordinate  £ and  i)  from 
x\  and  X2- 

The  vortcicity-vector  potential  method  was  validated  for  the  numerical 
simulation  of  180°  curved  duct  with  the  circular  section  and  the  transition 
simulation  in  a plane  channel  (Tokunaga  et  al  1991). 
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3.  Computational  Methods 


In  order  to  compute  accurately  the  shear  flow  turbulence,  the  fourth  order 
difference  method  is  used.  In  this  case,  the  first  and  second  derivative,  for 
example,  are  discretized  in  the  ^-direction  as 


du> 

d2u> 

w 


ijk 


(~^i-2,j,k  - 8w«+i,j,fc  - 8uji+i,j,k  + w»+2,j,fc)/(12  A £),  (5) 


ij,k 


— (— w»- 2,j,k  + ~ 30w,- J'k  + 16 Wj+ij-,*.  — 0Jj+2,j,fc) 


/(12AJ2), 


where  A£  represents  the  grid  interval  in  the  ^-direction  The  derivative  in 

rj-  and  ^-direction  as  well  as  the  cross  derivative  are  also  differenced  in  the 

same  way,  and  finally  the  set  of  the  ordinary  differential  equations  (ODEs) 
is  led. 

f = F(ui),  (6) 

<3  = (<^1, 2,1^2, 2,1,  (7) 

where  /,  J and  K denote  the  number  of  grid  points  in  £,  rj  and  z direction. 
For  time-stepping  in  the  set  of  ODEs,  we  apply  the  TVD  Runge-Kutta 
method  with  the  third  order  accuracy. 

The  same  discretization  procedure  is  adopted  for  the  Poisson  equation 
for  the  vector  potential,  where  the  point  Jacobi  is  used  in  £ and  z direction 
while  the  SOR  method  is  only  applied  in  tj  direction,  since  the  computation 
is  performed  by  Fujitsu  FACOM  VPP-800/63  of  the  architecture  of  pipeline 
with  the  huge  length.  The  convergence  of  the  Poisson  equation  gives  a sig- 
nificant influence  on  the  turbulent  field  and  insufficient  relaxation  weakens 
the  turbulence,  so  that  we  adopt  the  five  hundreds  relaxations  as  maximum. 


4.  Transition  Simulation  in  Channel  Flow 

In  order  to  validate  the  present  computational  method,  the  transition  pro- 
cess in  channel  flow  is  simulated.  The  Renolds  number  is  Re  = 8, 000,  and 
the  initial  condition  consists  of  the  Poiseuille  flow  with  2%  amplitude  2-D 
and  .02%  amplitude  3-D  T-S  waves,  which  is  the  standard  case  (Germano 
et  al  1991). 

The  time  evolution  of  the  wall  shear  stress  is  depicted  in  Fig.  2 (a), 
where  the  wall  shear  is  depicted  both  on  the  lower  and  upper  wall.  The 
transition  occurs  at  t = 180  due  to  the  impuritt  of  the  initial  condition. 
The  figure  2 (b)  shows  the  streamwise  velocity  distribution  average  in  the 
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Time 


(b)  Averaged  streamwise  velocity 
distribution 


(a)  Time  evolution  of  skin  friction 

@@ 


Figure  2.  Time  evolution  of  skin  friction  on  channel  wall  srafaces  (a)  and  streamwise 
velocity  distribution  averaged  in  spanwise  direction  (b) 


spanwise  direction,  which  shows  that  the  turbulent  flow  fully  develop es. 


5.  Numerical  Simulation  of  Channel  Flow  with  Two  Rows  of 
Transverse  Riblet 

Since  the  possiblity  of  drag  reduction  has  been  already  suggested  in  the 
single  row  of  transverse  riblet  in  the  periodic  domain  (Tokunaga,  1999), 
we  investigate  the  drag  reduction  in  two  rows  of  transverse  riblet  which  is 
shown  in  Fig.  3. 

The  Reynolds  number  Re  = Umax2S fv  is  chosen  as  8,000  and  the  height 
of  the  transverse  riblet  k = — 0.03  which  are  the  same  in  LES  of  the  single 
row  of  transverse  riblet. 


Figure  3.  Transverse  riblet 
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The  initial  turbulent  field  is  initiated  by  adopting  incidently  the  fully 
developed  channel  turbulent  flow,  which  is  obtained  in  the  previous  chapter, 
to  the  channel  with  two  transverse  riblets  in  order  to  elucidate  the  robust- 
ness of  the  present  computational  method.  The  computational  time  is  17 
hours  to  proceed  10  non-dimensional  time  by  making  use  of  single  vector 
processor  of  Fujitsu  VPP-800/63  which  is  the  parallel  super  computer  with 
63  processors  at  the  maximum  performance  504  GFlops  . Figure  4 (a)  and 
(b)  depict  the  time  history  of  the  total  drag  at  the  lower  and  upper  wall,  as 
well  as  the  skin  friction  and  pressure  drag  at  the  lower  wall.  In  the  present 
paper,  the  total  drag  denotes  the  sum  of  the  skin  friction  and  pressure  drag. 

It  is  shown  that  the  darg  on  the  flat  upper  wall  exceeds  the  total  drag 
on  the  lower  riblet  surface  after  the  initial  channel  turbulent  field  fits  to  the 
geometry  of  two  rows  of  transverse  riblet.  It  is  also  shown  that  the  form 
drag  is  too  small  compared  with  the  wall  shear.  From  this  point  of  view,  it 
is  concluded  that  the  transverse  riblet  has  potential  as  the  device  of  drag 
reduction. 
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(a)  total  drag  on  riblet  and  flat  wall 


Figure  4 • Time  history  of  total 


O 


Time 

(b)  skin  friction  and  form  drag  on 
riblet  wall 

(a),  skin  friction  and  form  drag  (b) 


6.  Turbulent  Field  on  Riblet  Surface 

In  this  section,  we  investigate  the  mechanism  of  drag  reduction  from  visual- 
izing the  turbulet  field  on  the  riblet  surface.  Figure  5 (a)  shows  the  contours 
of  the  spanwise  vorticity,  where  it  is  found  that  the  turbulet  vortices  are 
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attracted  to  the  transverse  riblet  and  weakened.  The  effect  is  the  same  as 
tiie  suction  used  in  the  drag  reduction  device. 

In  Fig.  5 (b),  we  show  the  contours  of  the  velocity  normal  to  the  wall 
which  is  averaged  in  the  spanwise  direction.  It  is  depicted  the  average  veloc- 
ity directs  to  the  riblet  surface.  Therefore,  it  is  thought  that  the  transverse 
riblet  has  the  same  effect  as  the  suction. 


(a)  Contours  of  spanwise  vorticity  (b)  Contours  of  the  normal  velocity 

@@ 

Figure  5.  Contours  of  spanwise  vorticity  on  x — y plane  (k  = 9 ) at  t = 110  (a)  and 
velocity  normal  to  the  wall  average  in  spanwise  direction  (b) 

The  iso-surface  of  the  spanwise  vorticty  is  depicted  in  Fig.  6.  It  is  found 
that  the  streak  structure  on  the  upper  flat  surface  is  smooth  and  strong, 
while  the  one  on  the  riblet  surface  is  blocked  on  the  two  rows  of  riblets  and 
weakened. 


Figure  6.  Iso  surfaces  of  streamwise  vorticity 
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7.  Conclusions 

The  complete  fourth  order  accurate  method  is  presented  for  LES  in  the 
vorticity-vector  potential  formulation  and  LES  is  performed  for  channel 
flow  with  two  lows  of  transverse  riblets.  As  results,  the  following  conclusions 
are  obtained. 

(1)  It  is  shown  that  the  channel  transition  simulation  is  successfully  carried 
out  by  the  present  computational  method. 

(2)  It  is  found  that  the  total  drag  is  reduced  after  the  turbulet  field 
developes  fully  on  the  riblet  surface. 


References 

Bobbit,  P.  J.,  Harvey,  W.  D.,  Harris,  C.  D.  and  Brooks,  Jr.,  C.  W.  (1991),  Natural 
Laminar  Flow  and  Laminar  Flow  Control,  eds.,  Barnwell,  R.  W.  and  Hussaini,  M., 
Y.,  Springer-Verlag,  , pp.  247-411 

Choi,  H.,  Moin,  P.  and  Kim  J.  (1994),  J.  Fluid  Mech.,  Vol.  262,  pp.  75-110 

Cortelezzi,  L.,  Lee,  K.  H.,  Kim,  J.  and  Speyer,  J.  L.  (1999),  Int.  J.  Comp.  Fluid  Dyn., 
Vol.  11,  pp.  79-92 

Germano,  M.,  Piomelli,  U.,  Moin,  P.  and  Cabot,  W.  H.  (1991)  Phys.  Fluids,  Vol.  29, 
pp.  1760-1765 

Liu,  Z.,  Liu,  Z.,  Liu,  C.  and  McCormick,  S.  (1994),  Int.  J.  Num.  Meth.  Fluids,  Vol.  19, 
pp.  23-40 

Mansour,  N.  N.,  Moin,  P.,  Reynolds,  W.  C.,  Ferziger,  J.  H.  (1979),  Turbulent  Shear  Flows 
/,  Springer-Verlag,  Berlin,  pp.  386-401 

Manuilovich,  S.  V.  (1994),  Theor.  Comput.  Fluid  Dynamics,  , pp.  31-47 

Mochizuki,  S.,  Kameda,  T.  and  Osaka,  H.  (1997),  Proc.  of  Int.  Conf.  Fluid  Engg.,  JSME 
Centenial  Grand  Conf.,  Tokyo,  JSME,  No.  97-203,  pp.  1655-1666 

Shimada,  M.,  Tokunaga,  H.,  Satofuka,  N.  and  Nishida,  H.  (1991),  JSME  Int.  J.  , Vol.  34, 
pp.  109-114 

Tokunaga,  H.,  Yoyeda,  K.  and  Satofuka,  N.  (2000),  Proc.  of  AIAA  10th  Computational 
Fluid  Dynamics  Conference,  Honolulu,  HI,  AIAA,  , pp.  937-946 

Tokunaga,  H.  (1999),  37th  AIAA  Aerospace  Sciences  Meeting  and  Exibit,  Reno,  NV, 
AIAA  Paper  99-0424,  , pp.  1-10 

Tokunaga,  H.  (1992)  Annual  Research  Briefs,  Center  for  Turbulence  Research,  Stanford 
Univ./NASA  Ames,  pp.  171-184 

Tokunaga,  H.  and  Yamauchi  (1997),  Int.  J.  Num.  Meth.  Fluids,  Vol.  25,  pp.  1107-1117 

Verzicco,  R.,  Yusof,  J.  M-.  Orlandi  P.  and  Haworth,  D.  (2000),  AIAA  J.,  Vol.  38,  pp.  427- 
433 

Weber,  C.  and  Ducros,  F.  (2000),  Int.  J.  Comp.  Fluid  Dyn.,  Vol.  13,  pp.  327-355 


